GS <- mutate(GS,Date = as.Date(Date, format = "%d-%b-%y"))
GSxts <- tk_xts(GS)
## Warning in tk_xts_.data.frame(data = data, select = select, date_var =
## date_var, : Non-numeric columns being dropped: Date
## Using column `Date` for date_var.
allDates = index(GSxts)
firstDate <- min(allDates)
lastDate <- max(allDates)-30 #find the last start_date
while(!lastDate %in% allDates)
  lastDate <- lastDate-1

result <- data.frame(`StartDate` = as.Date(character()), `OptionPnL` = double(), `HedgingPnL` = double(), `FinalPnL` = double(), `MaxDrawdown` = double(), `SharpeRatio` = double(),`StartPrice`= double(), `EndPrice` = double(), `AvgPrice` = double(), `avgGrowthRate` = double())

startD <- firstDate
for(startD in firstDate:lastDate){
  startD <- as.Date(startD)
  if(startD %in% allDates){
    endD <- startD+30
    #adjust the end date backwards if end date (a calendar day) is not in the xts
    while(!endD %in% allDates)
      endD <- endD-1
    
    xts_obj <- GSxts[paste(c(startD,endD),collapse = "/")]
    
    quantity = 100
    dates <- index(xts_obj)
    start_date <- min(dates)
    end_date <- max(dates)
    start_price <- as.numeric(xts_obj[start_date, "Close"])
    start_volatility <- as.numeric(xts_obj[start_date, "IV30"])
    
    df <- tibble(Date = dates)
    df$Close <- coredata(xts_obj[, "Close"])
    df$IV30 <- coredata(xts_obj[, "IV30"])
    avgChange <- as.numeric(mean(xts_obj[, "Change"],na.rm=TRUE))
    r <- 0.8 / 100
    X <- start_price/(exp(qnorm(0.25)*start_volatility/100*sqrt(30/365) - (r+0.5*(start_volatility/100)^2)*30/365))
    #sigma = start_volatility
    
    # Vary S and Time everyday
    #S <- df$Close
    #Time <- (end_date - df$Date) / 365
    
    #GBSOption(TypeFlag, S, X, Time, r, b, sigma)@price
    
    df_opt <- rowwise(df) %>%
    #this is the premium for one unit of call option  
    mutate(premium = GBSOption(TypeFlag = "c",
    S = Close,
    X = X,
    Time = as.numeric((end_date - Date) / 365),
    r = r, # interest rate
    b = 1.85/100, # dividend yield obtained from https://www.dividend.com/dividend-stocks/financial/investment-brokerage-national/gs-goldman-sachs/
    sigma = as.numeric(start_volatility/100))@price,
    
    #this is the delta of a call option (before negation)
    delta_hedge = GBSGreeks("delta", TypeFlag = "c", 
                            S = Close, 
                            X = X, 
                            Time = as.numeric((end_date - Date) / 365), 
                            r = r, 
                            b = 1.85/100, 
                            sigma = as.numeric(start_volatility/100))) %>%
    ungroup() %>%
      
    #delta hedging strategy selected: SHORT CALL LONG STOCK (from BlackS Scholes formula, such strategy should approximate a long position in risk free)
    mutate(Option_DoD_PnL = ifelse(Date == start_date, # On the 1st date, we count the cost of buying the option
    0, #quantity*premium, #on the first day, receive the call option premium and short the option
    -quantity*(premium - Lag(premium))), #if subsequently call option price rises, there is a loss
    
    Hedging_DoD_Pnl = ifelse(Date == start_date, 0, 
                             ifelse(Date == end_date, yes = quantity * Lag(delta_hedge) * (Close - Lag(Close)), 
                                    #at the last day, there is no rebalancing of number of shares.
                                    no = quantity * delta_hedge * (Close - Lag(Close)))), #long stock - if stock price increase, there is a profit
    
    DoD_PnL = Option_DoD_PnL + Hedging_DoD_Pnl) %>%
    mutate(PnL_to_date = cumsum(DoD_PnL),
           HPnL_to_date = cumsum(Hedging_DoD_Pnl), 
           OPnL_to_date = cumsum(Option_DoD_PnL))
    
    maxDrawDown <- {
    xs <- df_opt$PnL_to_date
    max(cummax(xs) - cummin(xs))
    } #the maxDrawDown could not be correctly calculated
    
    #The initial outflow of funds is the cost to buy stocks minus option premium received 
    #InitialInvt = (df_opt[[1,"delta_hedge"]]*df_opt[[1,"Close"]] - df_opt[[1,"premium"]])*quantity #OUTFLOW of funds
    #profitability = df_opt[df_opt$Date==end_date,"PnL_to_date"]/InitialInvt
    #df_opt<-mutate(df_opt, PortValue = InitialInvt + PnL_to_date, PortReturn = DoD_PnL/Lag(PortValue))
    
    #ggplotly(p=ggplot(df_opt) + geom_line(aes(TTM,Option_DoD_PnL),color = "blue") + ggtitle("option profit - TTM"))
    #ggplotly(p=ggplot(df_opt) + geom_line(aes(TTM,Hedging_DoD_Pnl))+ggtitle("stock profit - TTM"))
    
    #renderTable(tail(df_opt,1))
    
    
    #renderText(paste0("the Sharpe Ratio is ",round(SR,4)))
    #renderText(paste0("The maximum drawdown is ", round(maxDrawDown,4)))
    hedgingPnl <- as.numeric(df_opt[df_opt$Date==end_date,"HPnL_to_date"])
    finalPnl <- as.numeric(df_opt[df_opt$Date==end_date,"PnL_to_date"])
    optionPnl <- as.numeric(df_opt[df_opt$Date==end_date,"OPnL_to_date"])
    endPrice <- as.numeric(df_opt[df_opt$Date==end_date,"Close"])
    avgPrice <- as.numeric(mean(df_opt$Close,na.rm=TRUE))
    SR <- as.numeric((finalPnl/30)/stdev(df_opt$DoD_PnL, na.rm = TRUE)) #r is omitted and Sharpe ratio is calculated directly from the PnL
    result <- rbind(result,data.frame("StartDate" = start_date, "OptionPnL" = optionPnl, "HedgingPnL" = hedgingPnl, "FinalPnL" = finalPnl, "MaxDrawdown" = maxDrawDown, "SharpeRatio" = SR,"StartPrice"=start_price , "EndPrice" = endPrice, "AvgPrice" = avgPrice, "AvgGrowthRate" = avgChange))
    
  }}
    ggplotly(p = ggplot(GS) + geom_line(aes(Date, Close))) #stock close price
    ggplotly(p = ggplot(GS) + geom_density(aes(Close))) #density of close price
    ggplotly(p=ggplot(result) + geom_point(aes(AvgPrice,FinalPnL, color=AvgGrowthRate)) + ggtitle("avg price - final pnl"))
    ggplotly(p=ggplot(result) + geom_point(aes(AvgGrowthRate,FinalPnL))+ggtitle("avg growth rate - final pnl"))
    ggplotly(p=ggplot(result) + geom_point(aes(StartPrice,FinalPnL))+ggtitle("start price - final pnl"))
    ggplotly(p=ggplot(result) + geom_point(aes(EndPrice,FinalPnL))+ggtitle("end price - final pnl"))
    p1 <- ggplot(result) + geom_point(aes(AvgPrice, OptionPnL)) + ggtitle("avg price - option pnl")
    p2<-ggplot(result) + geom_point(aes(AvgPrice, HedgingPnL)) + ggtitle("avg price - hedging pnl")
    grid.arrange(p1,p2,nrow = 1)

    a1 <- ggplot(result) + geom_point(aes(AvgGrowthRate,OptionPnL))+ggtitle("avg growth rate - option pnl")
    a2 <- ggplot(result) + geom_point(aes(AvgGrowthRate,HedgingPnL))+ggtitle("avg growth rate - hedging pnl")
    grid.arrange(a1,a2, nrow = 1)

    a3 <- ggplot(result) + geom_point(aes(x = StartPrice, y = EndPrice, color = OptionPnL), alpha = 0.7)+ggtitle("start & end price - option pnl")
     a4 <-ggplot(result) + geom_point(aes(x = StartPrice, y = EndPrice, color = HedgingPnL), alpha = 0.7)+ggtitle("start & end price - hedging pnl")
     grid.arrange(a3,a4, nrow=1)

When we focus on hedging pnl, its range become larger when average price or average growth rate increases. When using 0.25 delta, we get a smaller final pnl, which is consistent with our expectation. It could also be observed that the option and hedging PnLs in a single trade with 25% delta fluctuate less violently or regularly against TTM as compared to ATM options in previous backtesting.

     kable(head(result,20))
StartDate OptionPnL HedgingPnL FinalPnL MaxDrawdown SharpeRatio StartPrice EndPrice AvgPrice AvgGrowthRate
2017-12-13 220.07269 128.730479 348.8032 348.8032 0.7580423 255.56 257.03 256.1119 -0.0309524
2017-12-14 210.73508 131.889557 342.6246 342.6246 0.7282751 255.48 257.03 256.1395 0.0735000
2017-12-15 204.73691 77.979243 282.7162 282.7162 0.6394765 257.17 257.03 256.1742 0.0815789
2017-12-18 227.06058 -6.320899 220.7397 220.7414 0.6540279 260.02 253.65 256.1125 -0.1760000
2017-12-19 232.23490 77.211706 309.4466 309.4467 0.7600284 256.48 250.97 255.6600 -0.4525000
2017-12-20 234.13694 125.935117 360.0721 360.0721 0.7916135 255.18 256.12 255.6420 -0.0180000
2017-12-21 232.49260 -71.733015 160.7596 160.7596 0.6032480 261.01 256.12 255.6663 0.0494737
2017-12-22 222.55065 -32.422465 190.1282 190.1282 0.5998324 258.97 256.12 255.3694 -0.2716667
2017-12-26 235.43117 156.499612 391.9308 391.9308 0.8726340 257.72 269.03 256.8571 0.4790476
2017-12-27 149.65258 534.855977 684.5086 684.5086 0.6720979 255.95 268.14 257.3533 0.4961905
2017-12-28 198.89344 460.331654 659.2251 659.2251 0.6348033 256.50 268.14 257.4235 0.6095000
2017-12-29 59.49473 651.501900 710.9966 710.9966 0.6556072 254.76 268.14 257.4721 0.6126316
2018-01-02 -308.80497 980.507094 671.7021 754.2604 0.5502747 255.67 272.23 259.9432 0.7940909
2018-01-03 221.18366 106.598993 327.7827 735.4081 0.1113464 253.29 260.04 260.1418 0.1986364
2018-01-04 223.80131 -106.387893 117.4134 813.7043 0.0236845 256.83 260.04 260.4681 0.3214286
2018-01-05 216.73966 -65.677868 151.0618 795.1294 0.0320384 255.52 260.04 260.6500 0.1605000
2018-01-08 226.99296 992.454914 1219.4479 1221.7959 0.4611731 251.81 257.10 260.1086 0.0718182
2018-01-09 235.04055 886.631393 1121.6719 1123.9790 0.4690417 253.94 246.35 259.8605 -0.2481818
2018-01-10 234.37942 873.655110 1108.0345 1108.0345 0.4917207 254.33 249.30 259.6495 -0.2109091
2018-01-11 235.13630 800.422280 1035.5586 1035.5586 0.4823469 255.13 249.30 259.9029 -0.2395238